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Abstract 

There has been devoted significant mathematical research to model the light prop- 
agation in tissue and to recover the absorption and scattering coefficients after pho- 
toacoustic inversion. Typically, the basic light propagation models considered there 
are the diffusion limits of the radiative heat transfer model. These equations are well 
suited for models where the elastic scattering is the dominant effect. If, however, the 
scattering is less pronounced, a single scattering approximation model is practicable. 
As we show in this paper, this approach is practically relevant in focused/sectional pho- 
toacoustic imaging. In this paper, we study analytical reconstruction formulas for the 
single scattering case. To realise the single scattering approach, we propose concrete 
physical experiments based on photoacoustical sectional imaging. 

Introduction 

In photoacoustic imaging (see e.g. [27, 17, 13, 24] for some mathematical and physical 
review papers), the interior of a small object is analysed by illuminating it with a short 
laser pulse and observing the acoustic wave which is hereby induced via the photoacoustic 
effect. The measurement of this pressure wave allows us (usually under the assumption that 
the acoustic wave travels with constant velocity and is free from any attenuation effects) to 
recover internal measurements in the form of the initially generated pressure 

p(°'(x)=7(x)MaWiW, i£E 3 , 

where 7 denotes the Griineisen parameter, which describes the change in pressure as energy 
is absorbed, /i a is the optical absorption coefficient of the material, and 3> is the light 
fluence. Often, the variations in 7 and § are neglected, so that P^ can be considered 
to be proportional to the absorption coefficient, which is used to characterise the material, 
see [9, 26] for some common reconstruction formulas for the initial pressure P(°\ 

Recent attempts have been made to model the light propagation in the material and 
to recover, so to speak in a second step, from the internal measurements P^ for different 
illuminations $ the absorption coefficient /i a and the scattering coefficient, which enters in 
the light propagation model, see e.g. [3]. The basic light propagation models considered 



*The work has been supported by the Austrian Science Fund (FWF) within the national research network 
Photoacoustic Imaging in Biology and Medicine, project S10505-N20. 



2 



Peter Elbau, Otmar Scherzer 



there are the diffusion limits of the radiative heat transfer model, the Boltzmann transport 
equation. These equations are well suited for models where the elastic scattering is the 
dominant effect. If, however, the scattering is less pronounced, a single scattering approx- 
imation model is practicable. This approach is practically relevant for focused/sectional 
photoacoustic imaging. In this paper, we study analytical reconstruction formulas for the 
single scattering case. Reconstruction formulas for the single scattering approach rely on a 
similar strategy as in the diffusion approximation [3, 4] and are based on deriving equations 
for quotients of independent measurement data. To realise the single scattering approach, 
we propose concrete physical experiments based on photoacoustical sectional imaging. 



1. Light Propagation Models 



Considering only elastic scattering, the light propagation inside the object can be modelled 
with a Boltzmann transport equation, the so-called radiative transfer equation, for the 
density of photons ip&(t, x) at the position x g R 3 at the time t € R moving in the direction 
d e S 2 of the form 

l d t ^{t,x) + {'dy x ^{t,x)) + l x t {x)^{t,x) = f"4^ f e«tf,0»v*(t,z)d*(tf), (i) 

c 4?r J S 2 

where the extinction or transport coefficent 

is given as the sum of the absorption coefficient /Lt a and the scattering coefficient fi s . More- 
over, c denotes the speed of light and O is the phase function, i.e. j- J Q Q({&, $)) ds($) is 
the probability that a photon heading into a direction 6 £1 C S 2 is scattered into the 
direction 

In photoacoustic imaging, however, the excitation happens with a short laser pulse with 
some fixed frequency v, and we are not interested in the exact light distribution as a function 
of time, but only in the total energy being absorbed at each point. So, let us switch in a 
first step to the energy fluence $^ originating from photons moving in the direction $ G S 2 
as new variable. We have the relation 

/oo 
ip&(t,x)cdt, x e R 3 , $ e S 2 , 
-oo 

where h denotes the Planck constant. 



Diffusion Approximation. In quantitative photoacoustic tomography, see e.g. [3] for a 

review, one commonly uses the diffusion approximation of this transport equation which 
takes the form 

div(aV§)(z) = A*a(a:)*(aO, *{x) = \ x € R 3 , (2) 

where 

Hx) = -^f $t(x)d8(0), i£E 3 , (3) 
4tt J S 2 

is the total light fluence at a point x and fi' s is the so-called reduced scattering coefficient. 
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To obtain this approximation, the dependency from the direction ■& is assumed to be at 
most linear. So we assume that the functions <E>,j and can for all x G R 3 and all tf,i)eS 2 
approximatively be written as 

S#(z)«fo(s) + <0,&(a:)) and 0«£,0» « Go + 0i<3,0) (4) 

for some (sufficiently smooth) functions 4>q : R 3 —¥ R, <pi : R 3 — > R 3 , and some constants 
O , ©l G R. Then we find for all x G R 3 and all <& G S* 2 that 

3 

(tf.V-^Cs)) pa (0,V&(a;)) + £ 0^*0*, = (0, Wofc)) + T D&(*)0 

and 

/ 0(^,^)^(2;) dsW«47r0 o ^o(a:) + ei / (0, t?) (0, 0i(sb)) ds(tf), 
,/S 2 JS 2 

where we used that the linear terms in d give zero when integrated over S 2 . To evaluate 
the integral therein, we switch to spherical coordinates with the polar angle a chosen as 
cos a = (•&,•&) and get for all £ S 2 , since the components orthogonal to •& vanish due to 
the axial symmetry of the integrand, that 

J ds(&) = J (d,tf) 2 dds(d) = f o cos 2 a sin a da^ = ytf. 

So, plugging the approximation (4) into the transport equation (1), which we average 
over time, we find for all x G R 3 that 

3 tr(D^i(x)) + ((ia,(x) + (is(x))<f>o(x) - ^ s (x)&Q(j)o(x) 

+ (•&, V0 o (x) + (n a { x ) + fi s (x))Mx) - y a e 1( f>i(x)) 

+ i9 T (D0 x (aj) - §tr(D& (x))) i9 « 0. (5) 

Ignoring the last term as it is of second order in ?? (since ■# T i9 = 1, we had to split off the 
trace term as it is of zeroth order) and using that 

©o«^-/ e((#,i?))ds(0) = 1 and fofx) a — f (x) ds(0) = x G R 3 , 

47T J S 2 47T „/ S 2 

we get by equating the coefficients in (5) the equation system 

|div^i(x) +/ia(x)$(x) w 0, ieR 3 ; 
V#(x) + (/i a (x) + ^4(x))^i(jb) « 0, z G R 3 , 

where ^(x) = (1 — |©i)/x s (x) is the reduced scattering coefficient in x £ R 3 . Plugging 
<f>i from the second equation in the first one, we arrive at the diffusion equation (2). A 
more general derivation for this diffusion equation, considering also higher order multipolc 
expansions of $^ and 0, can be found in [2]. 

Single Scattering Approximation. This diffusion model is well suited for materials 
where elastic scattering is the dominant effect. If, in contrast, the optical scattering is 
less pronounced, we may consider a single scattering model. In this case, we assume that 
scattered photons have a very low probability to interact again with the object, so that we 
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can simply treat them as lost and not track their movement anymore. We thus neglect the 
scattering term on the right hand side of the transport equation (1) and get for the light 
fluence <E>,j from photons moving in the direction ?9 £ S 2 the simplified transport equation 

(tf, V = -fit(x)^4x), x £ R 3 . (6) 

To specify boundary conditions for this differential equation, let us assume that no 
absorption or scattering occurs outside the object. Then the light fluence $,j is outside the 
object constant along all lines parallel to Thinking of the light fluence $^ being generated 
by a laser, which we for simplicity place infinitely far away in the direction — ■& £ S 2 , the 
fluence between the laser and the object should be known from the specifications of the 
laser. This means that we have for every i)eS 2 that 



lim <P»(y + \d) = <P { °\y) 7 y £ E*, (7) 

\—> — oo 

s(0) 



for some known inital light fluence $^ ' : E$ — > [0, oo) where 

E0 = {yeR 3 | (y,0) = O>. 

To formulate the solution of the differential equation (6) in a compact form, let us 
introduce for $ £ S 2 the notation 

x = xj} + x$d, x$ £ Ed, x-o £ R, (8) 

for the decomposition of a space coordinate x £ R 3 into a vector xjf £ R 3 orthogonal to 
$ and the component x§-d = (x, 1?) $ in direction of With the boundary data (7), the 
solution of the transport equation (6) is then given by 

(x) = $^ 0) (4 ) exp (4 + A<?) d\j 

for all x £ R 3 and all tfeS 2 . 

In particular, if we illuminate the object with just one laser beam of photons moving in 
a direction $0 £ S 2 , we have an initial light fluence of the form 

$$\x) = m(xl)6# (#), x£E#, (9) 

with the known initial total light fluence : E$ — s- R of the laser beam. Here, 5$ denotes 
the ^-distribution on the sphere at the point $0 £ S 2 defined by J g2 /(i9)<5^ (i?) ds(i9) = 
/(i?o) f° r a ll functions / £ C°°(S 2 ). In this case, the resulting total light fluence $, defined 
in equation (3), is given by 

§(x) = $^(4 )exp (- J *° Mt (4 + M )dX^j , x £ R 3 . (10) 



2. Reconstruction Formulas 



Let us now consider a photoacoustic measurement. I.e. we illuminate an unknown object 
resulting in a total light fluence <&. Then the photoacoustic effect generates an initial pres- 
sure 1 P^ ' which is proportional to the absorption coefficient /i a of the material and the 
total light fluence $ of the laser light: 

P {0 \x) = 7(i)/i a (x)$(x), ireR 3 . 

1 More precisely, we should refer to it as a pressure difference to the equilibrium pressure in the object. 
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The proportionality constant 7 is called the Griineisen parameter and describes the ther- 
modynamic properties of the material. 2 See e.g. [7] for a derivation of this relation. 

This initial pressure then initiates an acoustic wave propagating through the object. 
Usually, the simplified model of a homogeneous, elastic medium with constant speed of 
sound c s is made, which leads to the linear wave equation 

d tt P(t, x) = c 2 s A x P(t, x), t > 0, x € R 3 , 

d t P(0,x)=0, x e R 3 , (11) 

P(0,x) = P (0) (a;), xeE 3 , 

for the pressure P(t, x) at a point x g R 3 at time t > 0. 

There are a lot of articles discussing the problem of how to recover the initial pressure 
p(o) f rom 

some photoacoustic measurements m of the acoustic wave P outside the object. 
For classical measurements where the function m is of the form m : [0,oo) x dX — > R 
with m{t,x) = P{t,x) for some domain X containing the object, explicit reconstruction 
formulas could be obtained by Fourier methods [25, 28, 29] leading to the so-called universal 
back-projection formula [26, 19], which holds at least for the case where X is a half-plane, 
a cylinder, or an ellipsoid. Other approaches [9, 10, 15, 16, 14], reducing the problem to 
the inversion of the spherical means operator, also led to explicit reconstruction formulas 
for simple geometries X. And also from the more general setting of integral geometry used 
in [21, 22], reconstruction formulas for the photacoustic problem could be obtained. 

In the context of this paper, we simply want to assume that we are in some way able to 
recover from the photoacoustic measurements of the acoustic wave P outside the object the 
initial pressure P^ and thus have given the internal data 

P^\x)= 1 {x)^{x)^ i {x), xeR 3 , i = l,...,N, (12) 

for a certain number N G N of different illuminations of the object resulting in different 
total light fluences The aim would be to recover from these data the three material 
parameters: the Griineisen parameter 7, the absorption coefficient /i a , and the diffusion 
coefficient a (given by (2)) or the scattering coefficient fi s , depending on the choice of light 
scattering model. However, as it turns out, only two of these three parameters can be 
recovered as a function of the third from this sort of data. 



Diffusion Model. The focus of this work are explicit reconstruction methods for the single 
scattering approach. In comparison, the reconstruction for the parameters of the diffusion 
model (2) are discussed in the papers [1, 6, 3, 4]. To give an idea on what these methods 
are based, we shortly review the approach presented in [4]. 

The basis for explicit reconstructions there, as well as in our approach for the single scat- 
tering approximation, is to consider quotients of the known initial pressure data P^ , given 



by (12), for multiple illuminations i — 1, . . . , N. Indeed, from the diffusion equation (2), it 
jv 



follows (provided that P^\x) ^ for all x € R 3 ) that the quotients 



P$'(x) ®n(x) 



2 We have 7 = -zk-ac^ where a denotes the thermal expansion coefficient, c s is the speed of sound, and 
C p is the specific heat capacity at constant pressure. 
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fulfil for all x S R 3 and alH € {1, . . . , N — 1} the equation 

div(a$ 2 N V Ui )(x) = 0, i.e. ( ^$^ ,Vm(») ) = -A Ui (x). 



So, if at each point x € R 3 , the vectors Viti(x), i = 1, . . . , N — 1, span the whole space, we 
explicitly get the function V log(<7 < I > ^ r ) and can thus recover from some initial conditions (in 
principle from the knowledge of <J§ 2 N at one point) the function v given by 

v(x) = y/a(x)<t N (x), x e R 3 . 
From this function, we can then get the two combinations 

Pn = 7(>Wz) and Av{x) = fjjjz) Ay/a{x) ^ ^ 
v(x) y/a(x) v(x) a(x) yja{x) ' 

of the three material parameters 7, /i a , and er. However, as it is shown in [5], for given 
Dirichlet boundary data for all the physical parameters, these two combinations already 
uniquely determine the internal data P^ , so that no additional information can be extracted 
from the measurements. So, we can only express two of the three material properties as a 
function of the third, see again [5]. 

Single Scattering Model. For the reconstruction formulas in the single scattering model, 
it is enough to consider the initial pressures P^ and P^ generated by two different illu- 
minations to recover the material parameters. We choose to illuminate the object once with 
a laser beam from the direction —§ g S 2 and once with a laser beam from the opposite 
direction i3 e S 2 . According to formula (10) for the resulting total light fluences $1 and 
$2, respectively, we get with the notation (8) (using that = — x$ and x~t# = x$) for all 
x S R 3 the representations 

= i>^ 0) (:4)exp (- J \ t (xj + Ai?)dA^ and (13) 
$ 2 (x) = $( 0) (4)exp^-^ Mt(4 + Atf)dA^ (14) 

with some known initial total light fluences i^ -*, : E$ — > (0, 00) (for simplicity, we 
assume as in the diffusion model before that the whole space is illuminated). This we can 
plug into formula (12) to obtain the generated initial pressures 

p(°\x) = $ i °\xj) 1 {x)n a {x)exp(- J % t (4+Atf)dA^) and (15) 



P?'{x) = &f'(34)'r(x)n*(x)exp\- J nt(34 + \4)d\j, xeW. (16) 

Since only the combinations /it = Ma + Ms and 7/i a of the three unknown parameters 7, 
/i a , and /i s enter the expression for the initial pressure, there is no way to recover all three 
parameters, regardless of the number of measurements. Instead, we confine ourselves with 
reconstruction formulas for /i t and 7Ma- As in the diffusion model, the additional knowledge 
of any of the three parameters then immediately allows the reconstruction of the two others 
(at least in the domain where the initial pressures do not vanish). 
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Let us denote with Q the domain where the initial pressures do not vanish, i.e. 
n = {x e R 3 | P^\x)p( Q \x) > 0}. 
Then for x G fi, we find from (15) and (16) the identity 

$^(4)pf (x) _ 



log 



fit {xj + A#) dX — fi t (xj + X$) dX. 

Jx0 



^\xj)P^(x) 
Thus, taking the derivative in the direction we end up with 

1 



IM,{x) 



0, V K log 



for all x G fi. 



Assuming that 



(j, t (x)=0 for all xGR 3 \fi 



(17) 



(18) 



(meaning that whenever we have scattering or absorption at a point, this generates some 
pressure at that point), this completely determines the function fit- We can then get from 
this the product 7/z a with the relation 



7(x)^ a (z) 



\ ^ o) (4)*i o) 04) 



cxp Q J fj, t (xj + Ai9) dA^ , x G R 3 



(19) 



which follows directly from the formulas (15) and (16). 

However, we remark that measuring additionally the total light fluences 

$ ( ™\y) = lim + Ai?) and $^ oo) (y) = lim $ 2 (y - At?), y G (20) 

A— ^oo A— >-oo 

of the laser light behind the object, the exponential factor in the formula (19) can be 
recovered from these measurements by 



exp 



y G E$, i G {1,2}, 



which follows directly from (13) and (14). With this additional measurement, we do no longer 
need to know the function fi t at all points to calculate 7^t a (in contrast to formula (19)). 
In particular, we do not require the assumption (18) to determine the function jfj, a via the 
formula 



Pi°\x)P?\x) 



\ ^\xj)^\xj) \ <S>[°°\xj)^\xj) 



x G R' 5 . 



(21) 



3. Single Scattering Model in Photoacoustic Sectional Imaging 

A particular well suited example for the single scattering model is photoacoustic sectional 
imaging, see [11, 20] and [18, 23] for some experimental results with photoacoustic sectional 
imaging and [8, 12] for explicit reconstruction formulas for the initial pressure. 
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In this sectional setup, the object is not uniformly illuminated (as it is typically the case 
for standard photoacoustic imaging), but the laser light is focused so that only one slice of 
the object is illuminated. For simplicity, we want to assume that the illumination of the 
object is done with a single laser beam coming from a laser placed at infinity. As before, we 
consider the illumination from two opposite directions. And to further simplify the notation, 
we choose these directions to be — d = (—1, 0, 0) and $ = (1, 0, 0). Moreover, this beam shall 
be perfectly focused onto the illumination plane {(x,y, z) g R 3 | z — 0}, so that the two 
corresponding initial total light fluences $^°' ) and $2°' °f the laser beam, appearing in the 
equations (13) and (14), have the form 

$?\o,y,z) = ^ 0) (y)S(z), y,z e R, i G {1,2}, 

for some known functions $ 2 °' ) : R —> (0, oo). 

In practice, of course, some scattering effects will still occur, leading to a larger illu- 
mination region inside the object. To diminish this effect, people started to use focusing 
detectors for the measurement of the acoustic waves, see e.g. [11, 20]. These detectors are 
tuned in such a way that pressure waves originating from points outside of the desired illu- 
mination plane interfere destructively on the detector surface so that these waves contribute 
considerably less to the measurements than those originating from the illumination plane. 

From the modelling point of view, these focusing detectors simply suppress the detection 
of those pressure waves which are generated by the absorption of scattered photons (unless 
they are only scattered inside the illumination plane or multiple times in a way that they 
end up being absorbed in the illumination plane again). On the other hand, this absorption 
of scattered photons is exactly the effect which we neglect in the single scattering model. 
Therefore, the single scattering model seems to be a good approximation for the modelling 
of photoacoustic sectional imaging with focusing detectors. 

Remarking that for the direction i9 = (1,0,0), the decomposition (8) simply reads 
(x,y,z) = (0,y,z) + a;(l,0,0), we get from the equations (15) and (16) that the result- 
ing initial pressures P^ and in this single scattering model are given by 

pf\x,y,z) = Pf\x,y)5{z), x,y,z€R, ie{l,2}, 

where the functions P^ , P^ : R 2 — > R are defined by 

Pi Q) (x,y) = <b { °\y) 1 (x,y,0)fi a (x,y,0)exp(- J Mt(A, y, 0) dX^j and (22) 
p( Q) (x,y) = <b { ° ) (yh(x,y,0)fi a (x,y,0)exp(- J ^(X,y,0)dx\, (x,y)eR 2 . (23) 

Modelling the propagation of the pressure wave with the linear wave equation (11), 
we derived in [8] explicit reconstruction formulas for the initial pressures Pj; \ i <E {1,2}, 
in the illumination plane for different detector setups, where the focusing detectors are 
approximated by either point, integrating line, or integrating plane detectors. 

As the only difference between the expressions (22), (23) and (15), (16) is that the first 
ones do not depend on the distance z to the illumination plane, we can use the exact same 
derivation as we used to get the formulas (17), (19), and (21) to recover fi t and the product 
7/i a . We therefore find for fj, t that 

i P^fr v) 

IMt(x,y,0) = -d x log „ 2 (0) l ,y > for all fay) € SI, 



pr(x, y ) 
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where Cl = {(x,y) € R 2 | P^\x,y)P^\x,y) > 0}. And if we assume again that the 
transport coefficient /i t fulfils that fa(x,y, 0) = if (x,y) g R 2 \ fl, then this allows us to 
calculate 7/i a via 



j(x, y, 0)fi a (x, y, 0) = 



\| 6XP 



i(0)/ 



/xt(A,y,0)dA , (i,j)eK 2 . 



As in the previous section, this can also be simplified to 
7(x,y,0)L( a (x,y,0) 



(0), 



with the additional measurement of the total light fluences 

$ ( i °°\o,y,z) = $ ( i °°\y)6(z), y, z € R, i E {1 , 2}, 
behind the object as defined in (20) with ■§ = (1, 0, 0). 



(£,2/)€R 2 , 



Conclusion 



We have shown explicit reconstruction formulas for photoacoustic imaging for the three main 
physical parameters, the Griineisen parameter, the absorption and the scattering coefficient, 
in a single scattering light propagation model. Here, in analogy, what has been pointed out 
earlier [3, 4] for the diffusion model, it is also only possible to recover two of them as a 
function of the third. 

Moreover, we have argued that this single scattering model is a good approximation for 
photoacoustic sectional imaging where focusing detectors are used to measure only acoustic 
signals originating from the illuminated plane. 
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